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We use replica exchange Monte-Carlo simulations to measure the equilibrium equation of state of 
the disordered fluid state for a binary hard sphere mixture up to very large densities where standard 
Monte-Carlo simulations do not easily reach thermal equilibrium. For the moderate system sizes we 
use (up to N = 100), we find no sign of a pressure discontinuity near the location of dynamic glass 
singularities extrapolated using either algebraic or simple exponential divergences, suggesting they 
do not correspond to genuine thermodynamic glass transitions. Several scenarios are proposed for 
the fate of the fluid state in the thermodynamic limit. 



I. INTRODUCTION 

Simple liquids, crystals, glasses, powders, and colloidal 
dispersions are frequently modeled using hard spheres. 
Although considered as one of the simplest models in 
condensed matter physics, hard spheres exhibit a compli- 
cated phase behavior that is not fully elucidated. In par- 
ticular, a well-established first order fluid-solid transition 
exists in three dimensions. For monodisperse systems, it 
occurs from volume fractions ipi = irpa 3 /6 « 0.492 to 
ip s rs 0.545 Jp is the number density and a the particle 
diameter) [lHfl]. However, the metastable fluid branch 
persists for (p > <pf , and its fate at large <p remains a de- 
bated subject @. Since the fluid cannot exist above the 
maximum density of the cubic centered crystal structure, 
ipf cc w 0.74, it may either go unstable, or it may exhibit 
a singularity with a diverging (dimensionless) pressure, 
Z = f3P/p, and a vanishing (also dimensionless) isother- 
mal compressibility, x = 5p/6(f)P). Here (3 = (fc^T) -1 , 
with T the temperature and ks the Boltzmann con- 
stant. Additionally, a thermodynam ic g lass transition 
could possibly occur along the way , characterized 

by a diverging timescale for structural relaxation [l2l - [l5l ] , 
a change of slope in the equilibrium equation of state 
Z(ip), and a jump in the compressibility. These features 
would be the analog, for hard spheres, of the glass tran- 
sition observed in glassforming liquids, characterized in 
particular by a diverging viscosity and a jump in the spe- 
cific heat [16| . It is the aim of this work to search for a 
thermodynamic signature of the glass transition in hard 
spheres. 

Studying the metastable fluid branch by simulations 
is complicated since the system naturally tends to form 
the crystal phase [l7|, [HI, & t least in three dimensions 
[l9T | . Since pressure is very dependent on the existence of 
small amounts of crystal nuclei, excluding ordered config- 
urations from the sampling is critical to obtain the real 
pressure-density relationship [13, [Uj]- An efficient way 
to overcome the somewhat arbitrary exclusion of crys- 
talline states from the sampling is to introduce size poly- 



dispersity to avoid, or at least considerably delay, crystal 
formation. One must then work between several con- 
straints: polydispersity must be large enough to prevent 
ordering, but small enough that a qualitatively different 
physics, specific to very polydisperse systems, does not 
set in. For instance, phase separation can occur in mix- 
tures [20], or fractionation in systems with continuous 
polydispersity [2lj . These phenomena have counterparts 
even for disordered states, since multiple glass transitions 
might occur in polydisperse systems, where for instance 
large particles are arrested in a sea of small ones that 
still easily diffuse [22j]. In this work we use a 50:50 binary 
mixture of hard spheres with a diameter ratio 1.4, large 
enough to efficiently prevent crystallization, but which 
shows no sign of multiple glass transitions. 

The final problem to be overcome is also the most dif- 
ficult one: approaching the glass transition at thermal 
equilibrium is hard in systems where the viscosity be- 
comes large because the timescale to reach equilibrium is 
simultaneously diverging. On this aspect, numerical sim- 
ulations could potentially outperform experimental work 
since it is possible, at least in principle, to imagine algo- 
rithms that have no 'physical' counterpart but still allow 
a proper exploration of the configuration space, and thus 
of the thermodynamic properties of the system |23j . Sev- 
eral such 'smart' algorithms exist in various context of 
statistical mechanics, such as umbrella sampling which 
makes use of biased statistical weights, replica exchange 
or parallel tempering where copies of the system at var- 
ious thermodynamic states are run in parallel to avoid 
being trapped in free energy minima [24J-126J, or cluster 
and swap algorithms which implement unphysical parti- 
cle moves to speed up equilibration [2?], HI] . 

Although commonly used and very successful in many 
areas of condensed matter, such methods have compar- 
atively been much less used in numerical studies of the 
glass transition, for several reasons. Firstly, the glass 
transition is mostly defined by, and studied via, dynamic 
properties, and so it is vital to use physical microscopic 
dynamics, which inevitably yields slow dynamics. There 
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are nevertheless interesting thermodynamic properties to 
be investigated in glassforming materials, for which im- 
plementation of particle swaps [Hj], cluster moves [30l ]. 
Wang-Landau sampling [3ll l32j | , or parallel tempering 
[33M36l | have all been implemented. Of course these differ- 
ent methods can be combined to improve further the effi- 
ciency. This has led in particular to strong claims about 
both the absence [3(| and presence [2(| [37[ of thermo- 
dynamic glass transitions in various glassy fluid models 
(including hard spheres), but also raised debates about 
the real efficiency of the various numerical alg orithms to 
study systems with slow dynamics [33j-[s^ [38|. 

Here, we employ the replica exchange Monte-Carlo 
(REMC) method jH H, |39|, HI, as recently adapted 
to systems composed of hard particles [|| . The idea is to 
simulate several replicas of the same system at different 
but close enough thermodynamic states to allow efficient 
exchanges between the replicas 23]. For soft interpar- 
ticle potentials, the most common ensemble expansion 
is that performed in temperature where each replica fol- 
lows a canonical ensemble simulation and the ensembles 
are set at different temperatures. To take advantage of 
the REMC algorithm for hard spheres, one needs to ex- 
pand the isobaric-isothermal ensemble in pressure [4lj |. 
each replica evolving at a different pressure @. Repli- 
cas having the larger pressures can escape from locally 
stable free energy minima through successive exchanges 
with replicas at lower pressures 42j. 



Using REMC, we have been able to reach thermal equi- 
librium for hard spheres up to very large densities where 
standard Monte-Carlo algorithms do not allow proper 
sampling of the configuration space [H, [44[ . We were 
thus able to study the thermodynamic properties of the 
disordered fluid branch of a binary hard sphere mixture 
over a broad density range which includes both the mode- 
coupling, ip mc t, and Vogel-Fulcher-Tamman, ip v [ t > v?mct, 
dynamic singularities finding no thermodynamic signa- 
ture for any of them, at least for the moderate system 
sizes we used, up to N = 100. While the absence of a 
genuine transition at ipmct can be established by standard 
numerical methods [H, HH , it is the main new result 
of this work that the same phenomenon seems to occur 
also at ip v ft- 

The paper is organized as follows. In Sec. [TT] we 
describe in more detail the model we use, and review 
the various 'critical' volume fractions that have been re- 
ported in previous work. In Sec. Mil we provide details 
about the REMC simulations. In Sec. IIVI we perform 
several tests to ensure that a proper sampling of config- 
uration space has been done. In Sec. [V] we describe our 
equilibrium results for the thermodynamics of the sys- 
tem. In Sec. I VII we investigate even higher densities, for 
which thermal equilibration could not be reached. Fi- 
nally, we discuss our results in Sec. IVHI 



Definition 


Volume fraction 


Onset of glassy dynamics 
Mode-coupling theory, Eq. |T} 
Vogel-Fulcher-Tamman, Eq. (J2]| 
Dynamic scaling, Eq. ([3| 
Diverging pressure (lower bound) 


Vonsct ~ 0.56 
fmct = 0.592 
ipvft = 0.615 
ipo = 0.635 
V?iow = 0.662 



TABLE I: Values of the relevant volume fractions character- 
izing the physical behaviour of the fluid for the binary hard 
sphere mixture studied in this work. 



II. CRITICAL DENSITIES IN A BINARY HARD 
SPHERE MIXTURE MODEL 

Previous work on binary mixtures suggests that a 50:50 
binary mixture of hard spheres with a diameter ratio of 
1.4 is a very efficient way to prevent crystalline ordering 
even at large densities 0, 0, E^j. We will use N — 
Na + Nb particles, Na and Nb denoting the number of 
small and large particles in the mixture, respectively. We 
work in units where the diameter of the small particles 
is unity, a aa = f • 

Moreover, the dynamics of small and large particles is 
strongly coupled so that the slow relaxation and location 
of the putative dynamic glass singularities yields consis- 
tent results for both components of the mixture [H, |46| . 
Thus, this model seems well-suited for investigating the 
existence of a thermodynamic glass transition of hard 
spheres. In a previous (bidimensional) study where effi- 
cient cluster Monte-Carlo moves were used (30j |. a very 
large polydispersity was introduced, with the unwanted 
result that large particles seemed to arrest at a density 
where small particles could still easily diffuse, making the 
identification of dynamic singularities somewhat ambigu- 
ous [1|. 

Previous numerical explorations of the dynamics of the 
present binary mixture revealed the existence of very slow 
dynamics and possible dynamic singularities at large vol- 
ume fraction, with no interference from the crystalline 
phase [lj| . Several relevant values of the packing frac- 
tions have been reported using different definitions and 
theoretical approaches, and we summarize them in Table 
III 

First, the dynamics of the system slows down and 
starts to become non-exponential above (p onS ct ~ 0.56, 
which can thus be seen as the onset density for slow dy- 
namics in this system. 

Second, the location of several dynamic 'singularities' 
can be defined and have been numerically studied. An 
algebraic divergence of the relaxation time, 



(^rnct - (p)~ 



(1) 



as predicted by mode-coupling theory [47[ , can be located 
near ip mc t ~ 0.592. However, simulations also revealed 
this density to be a crossover since the equilibrium re- 
laxation time can be measured at and above (/'met where 
it remains finite [HI, [46j] . This suggests that a different 
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functional form should be used to extrapolate a possible 
divergence of the relaxation time. 

A popular functional form for r(ip) is the so-called 
Vogel-Fulcher-Tamman (VFT) expression (l6j . 

r ~ ^ exp ( — ] , (2) 

V^vft - <PJ 

which yields, for the present system, the value ip v ft~ 
0.615, A and t m being additional fitting parameters }15| . 
As opposed to the mode-coupling singularity, standard 
simulations fail to access such a large packing fraction 
in equilibrium conditions, since the largest state point 
investigated in Ref. [l5[ is tp = 0.597 < <p v ft- 

Using a combination of scaling arguments involving 
both direct simulations of hard particles, and a soft har- 
monic repulsion at very low temperatures, recent numer- 
ical work provided support for the existence of a slightly 
different, stronger dynamic divergence EE BEE El, 

with the preferred values S « 2.2 and <po w 0.635. 

Finally, taking the view that no thermodynamic glass 
transition occurs, one must conclude that dynamics 
should arrest when particles come into contact and no 
particle move can take place. In this perspective, r must 
diverge simultaneously with the pressure Z at the ran- 
dom close packing or jamming density, V'rcp, which can 
then be empirically defined as the end point of the equi- 
librium equation of state of the fluid branch [4i| . In prac- 
tice this is hard to measure because the system falls out of 
equilibrium and becomes a nonergodic hard sphere glass 
much before getting to jamming, such that only lower 
bounds to the location of the diverging pressure can be 
numerically determined. For the present system, previ- 
ous work reported the value </?i ow ~ 0.662 as the tightest 
lower bound on </J rC p, obtained by rapid compressions of 
carefully equilibrated fluid states [3, This result 

indicates that the putative end point of the metastable 
fluid branch for this system is above tp\ ow — 0.662. 

III. THE REPLICA EXCHANGE 
MONTE-CARLO METHOD 

The partition function in the extended ensemble stud- 
ied in the replica exchange Monte-Carlo method we use 
is given by 0, l4lj 

n r 

Qcxtcndcd = Y\_ QnTP, , (4) 

where QntPi is the partition function of the isobaric- 
isothermal ensemble of the system at pressure Pi, tem- 
perature T, particle number N. The important new pa- 
rameter is n r , the considered number of replicas of the 
system. 



This extended ensemble is sampled by combining stan- 
dard NTPi simulations on each replica (involving both 
trial displacements of single particles and trial volume 
changes) and replica exchanges (swap moves at the 
replica level). To satisfy detailed balance, these swap 
moves are performed by setting equal all a priory prob- 
abilities for choosing adjacent pairs of replicas and using 
the following acceptance probability [g, .41] 

P acc = min(l,exp[/3(P - P,-)^ - Vj)]), (5) 

where Vt — Vj is the volume difference between replicas 

1 and j. Adjacent pressures should be close enough to 
provide nonnegligible exchange acceptance rates between 
neighboring ensembles. In order to take good advantage 
of the method, the ensemble at the smaller pressure must 
also ensure large jumps in configuration space, so that the 
larger pressure ensembles can be efficiently sampled. 

The probability for selecting a particle displacement 
trial, Pd, for selecting a volume change trial, P Vl and a 
swap trial, P s , are fixed to 

P d = n r N/(n r (N + l) + w), 

P v = n r /{n r {N + l) + w), (6) 

P s = w/(n r {N + 1) + w), 

where w -C 1 is a weight factor. Note that Pd + P v + P s = 
1, as it should. The probability density function to have 
the next swap trial move at the trial n t is given by 

P(n t ) = P s exp(-P s n t ). (7) 

Hence, one may obtain the next swap trial move from 
n t = — ln(£)/P s , with £ being a random number uni- 
formly distributed in the interval ]0, 1] [U [52j. We set 
all particles of a given replica to have the same a priori 
probability of being selected to perform a displacement 
trial. The same is true for selecting a replica for perform- 
ing a volume change trial. 

The trials [l,n t — 1] are displacements and volume 
changes, and so, they can be independently performed on 
the replicas. This has the advantage of being easily par- 
allelized. The algorithm is parallelized in four threads, 
since quad core desktops are used, but could be more ef- 
ficiently parallelized in n r threads. Since all swap trials 
are performed in a single core, the efficiency of the par- 
allclization increases with decreasing w. We employed 
w = 1/100. Verlet lists are used for saving CPU time, 
which can be quite large for the replicas evolving with 
the highest pressure values. 

Our simulations are performed in two steps. All simu- 
lations are started by randomly placing particles (avoid- 
ing overlaps), so that the initial volume fraction is cp = 
0.30. We first perform about 2 x 10 13 trial moves at the 
desired state points, during which we observe that the 
replicas reach a stationary state. We then perform more 

2 x 10 13 additional trials during which various measure- 
ments are performed, with results described in the fol- 
lowing sections. 




q> q> 

FIG. 1: a) qaa{(Taa) (squares), qab{oab) (diamonds), and 
SflB(crss) (circles) as a function of (p. The total pressure, 
Eq. (O, is shown as triangles. It agrees well with the set pres- 
sure values (light bullets), b) The number of AA (squares), 
AB (diamonds), and BB (circles) neighbors as a function of 
ifi with bullets indicating the total number of neighbors per 
particle. All data correspond to N — 100. 



The maximum particle displacements and volume 
changes for trial moves are adapted for each pressure to 
yield acceptance rates close to 0.3. Thus, particle dis- 
placements and volume changes of ensembles having high 
pressures are smaller than those associated to ensembles 
having low pressures. An optimal allocation of the repli- 
cas should lead to a constant swap acceptance rate for 
all pairs of adjacent ensembles. For a temperature ex- 
pansion, the efficiency of the method peaks at swap ac- 
ceptance rates close to 20% (53] . In this work, we use 
instead a geometric progression of the pressure with the 
replica index. In Sec. |V] we report results for various 
system sizes, N = 60, 80, and 100 using n r — 14, with 
fiP varying from 38 to approximately 5.8, the geometri- 
cal factor being 0.865. In Sec. |VI] we present additional 
results where the largest pressure is f3P = 100, N = 60, 
n r = 18, and the geometrical factor is 0.840. 



IV. THERMALIZATION TESTS 

The aim of this work is to provide new, reliable thermo- 
dynamic information at large densities where thermaliza- 
tion becomes a severe issue for standard algorithms. This 
means in particular that the algorithm must be able to 
sample accurately a phase space where ergodicity is po- 
tentially broken in the thermodynamic limit. Such severe 
sampling conditions are also met for instance in systems 
such as spin glasses 40] . It is crucial to establish whether 
the produced results are indeed representative of thermal 
equilibrium, as we now discuss. 

As a first check we verify that the pressure measured 
from the configurations sampled by the replicas in each 
ensemble yield results consistent with the values set nu- 




trials 



FIG. 2: Random walk in pressure space for N = 100 and 
n r = 14 for three arbitrarily chosen replicas. Thermodynamic 
state '1' corresponds to the highest pressure and '14' to the 
lowest. All replica visit several times both lowest and highest 
pressure states. 

merically. Pressure and structure are related by [20| 

j3P 2np v v o . . 

p ^ 

a 7 

where a and 7 run over species A and B, and x a , <r Q7 , 
and g aj respectively being the fraction of particles in 
species a, the contact distance between a and 7, and 
the partial radial distribution functions of species a and 
7. Note that <7 a7 (<7 Q7 ) must be evaluated using a careful 
extrapolation of g a ~t(r) towards contact. Thus, we may 
split the excess pressure into three contributions, corre- 
sponding to the AA, AB, and BB interactions. These 
contributions are shown in Fig. [TJa together with the to- 
tal pressure obtained from Eq. ©. As can be seen, the 
measured pressure agrees very well with the values im- 
posed numerically. Furthermore, a smooth behavior is 
obtained for all g a .y(a ai ) as a function of ip suggesting 
that adequate sampling has been performed. 

For all ip the largest contribution to the excess pres- 
sure is that of the large-large pairs (BB), followed by 
the large-small (AB) and the small-small (AA) pairs, in 
that order. In the right panel of Fig. [1] the evolution of 
the average number of neighbors is shown, obtained by 
integration of the partial pair correlation functions in a 
spherical shell of constant thickness 0.2<jaa- The total 
number of neighbors per particle is also shown. The num- 
bers of neighbors are consistent with the contributions to 
the excess pressure, i. e., they increase following the or- 
der AA, AB, and BB at all (p and they increasing with 
ip, albeit more slowly than the pressure. This satura- 
tion is physically expected since the number of neighbors 
remains finite even when the pressure diverges near jam- 
ming. 

We mentioned in the introduction that preventing crys- 
tallization is in principle dealt with by using a binary 
mixture. Since evidence for this stems from standard nu- 
merical approaches, it remains to be seen whether REMC 
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FIG. 3: Probability distribution functions (PDFs) of volume 
fraction fluctuations for each of the n r = 14 pressure val- 
ues, N = 100. All distributions are stationary, featureless, 
and symmetric and have sufficient overlap to allow for replica 
exchanges. 



finds more easily the crystalline phase or not. Indeed pre- 
vious work on the monodisperse system showed that the 
REMC algorithm is not only capable of forming the crys- 
tal phase but also to accurately predict the liquid-solid 
transition Thus, if a crystal is the preferred state, 
we expect to see signs of local orientational order. We 
checked this by computing the well-known order param- 
eter Qq, as denned for instance in Refs. [HEH, [54|, which 
is very sensitive to any trace of local angular order [l8[ . 
We evaluate it separately for AA, BB, and AB pairs. In 
all cases and at all densities, Q$ is very close to the value 
of a completely random system of points 18]. Moreover, 
the three Qe values do not evolve significantly during 
the runs. Thus, we can safely conclude that if the crystal 
phase corresponds to the equilibrium state of this partic- 
ular binary mixture, it is sufficiently metastable not to 
affect our results regarding the disordered state. 

For replica exchange methods to provide an efficient 
sampling of phase space, it is important to check whether 
all simulated replicas visit the entire set of thermody- 
namic conditions several times. This is a necessary con- 
dition for thcrmalization because this ensures that the 
configurations contributing to the thermodynamic aver- 
ages are very different as the low pressure replicas evolve 
rapidly and have large displacements in configuration 
space. In Fig. [5] we show the evolution of three randomly 
selected replicas making a random walk among the dif- 
ferent pressure states. We observe that all replicas con- 
tribute several times during the course of the production 
run to both the highest and the lowest pressure states. 
With the imposed geometric progression of the pressure, 
we find that the acceptance rate has a small drop near 
ip w 0.58, which suggests that thcrmalization becomes 
much harder above these densities. It is intriguing that 
this corresponds roughly to tp mc t, above which it also be- 
comes hard to reach thermal equilibrium using standard 
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FIG. 4: a) Equation of state, Z = /3P/p = Z(ip), for N = 60 
and increasing number of replicas from n r — 1 to n r — 14. 
Only data for n r > 8 are reproducible, while data for a 
smaller n r are not thermalized. b) Checking the FDT re- 
lation, Eq. ([9]), for density fluctuations. The lines show 1/x 
obtained from taking the derivative of the equation of state, 
while symbols are direct evaluation using spontaneous fluctua- 
tions of the density. The FDT holds with good accuracy both 
for equilibrated systems (n r = 14, bottom) and for nearly 
frozen ones (n r = 6, top). 



numerical tools. This numerical bottleneck suggests a 
faster decrease of the number of accessible configurations 
as pressure increases, or a faster increase of the barriers 
separating long-lived metastable states. 

Each replica evolves at a given predefined pressure 
Pi. Therefore, the volume Vi of replica j is a fluctuat- 
ing quantity, and it is interesting to focus on the prob- 
ability density functions (PDFs) of the volumes Vi, or 
equivalently of the corresponding volume fractions ifi. 
In cases where the system remains trapped in long-lived 
metastable states, the PDFs may be distorted or may 
contain peaks or shoulders which help detecting a lack of 
thermalization. Additionally, these features of the PDFs 
typically disappear as time increases and thus help re- 
vealing whether measurements are performed in station- 
ary states. We observe that the PDFs evolve in the first 
simulation steps but they then become both symmetric 
and take a Gaussian shape. The resulting functions are 
shown in Fig. [3] for N = 100 with low volume fractions 
PDFs correspond to low pressures. Notice that the PDFs 
have a larger peak and become narrower as pressure in- 
creases, which reflects the fact that the compressibility 
decreases. 

Ergodicity implies that the same results should be ob- 
tained independently of the set initial conditions and of 
the parameters of the simulation. We checked the re- 
producibility of our results by running simulations with 
n r = 1, 2, 4, 6, 8, 10, 12, and 14, for N = 60, for the 
same highest pressure (/3P = 38) and geometrical fac- 
tor of 0.865. Three independent runs were carried out 
with n r = 1, having all different initial conditions. For a 
fair comparison, all simulations lasted four weeks running 
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on identical single cores and no parallelization was imple- 
mented for this particular test. From the measured PDFs 
at each pressure, we measure the averaged density to cal- 
culate the equation of state Z(ip), which are reported in 
Fig. [4] We observe that results become reproducible only 
when n r > 8, which corresponds to simulations where 
the lowest pressure is below ip — 0.58 and yields ther- 
malized results. For smaller n r , Z is always larger than 
that obtained for n r > 8, suggesting that thermal equi- 
librium had not been reached. In particular, the three 
independent runs with n r — 1 are well above the equili- 
brated curve and distinct from one another. This implies 
that runs with n r < 8 are nonergodic and do not sample 
the configuration space accurately at large densities even 
with a large number of trials. In particular, this means 
that a standard Monte-Carlo algorithm would not yield 
equilibrium results at large density, and that it is clearly 
outperformed by the REMC simulation scheme we use in 
this work. 

A final test for equilibrium was suggested by Santen 
and Krauth (30j . At thermal equilibrium the sponta- 
neous fluctuations of density are related to the isother- 
mal compressibility, which is defined from the response 
of the pressure to an infinitesimal change in density in 
the linear response regime. This relation is thus a form 
of fluctuation-dissipation theorem (FDT), which is de- 
rived using the hypothesis that states are sampled with 
the equilibrium Gibbs measure: 



X = N 



(p 2 ) ~ (p? 
(p) 2 



dp 



5(J3P) ' 



(9) 



Checking whether this relation is satisfied by the data 
is therefore in principle a good way to check thermali- 
sation. In practice, this means checking the existence of 
a quantitative relationship between the broadness of the 
PDFs in Fig. [3] and the location of their averages. 

We followed these two routes for obtaining \. The 
fluctuations are directly measured from the simulations, 
while the response function is obtained by first fitting the 
pressure locally to a smooth polynomial function before 
taking the derivative with respect to density. In Fig. @] 
we present our results showing 1/x as a function of den- 
sity using both methods. When n r = 14 and results are 
reproducible, we find that the FDT relative to fluctua- 
tion density is well satisfied, which comes as an additional 
proof that our data are representative of thermal equilib- 
rium. However we note that for runs with a small num- 
ber of replica all concentrated in the high density regime 
which appeared far from equilibrium in the left panel 
of Fig. [H the FDT is also satisfied with a good accu- 
racy. In that case, all replicas belong to the glassy state 
and are basically frozen in a single 'basin' where they 
sample quasi-equilibrium short-lived fluctuations. Stud- 
ies of FDT violations in aging glasses have indeed shown 
that deviations from FDT appear only when considering 
those deg rees of freedom that relax very slowly in the 
glass 55]. This suggests that the FDT test suggested 
in Ref. [30j is only effective in a narrow density regime 



where a complete separation of timescale does not make 
Eq. @ valid even very far from equilibrium. That is, the 
test seems useful to detect a slow evolution and thus, the 
consistency of both \ determinations only guarantees a 
stationary state has been reached. This is a necessary 
but insufficient condition for equilibrium. 



V. THERMODYNAMIC RESULTS AT 
EQUILIBRIUM 

In previous sections, we provided evidence that the 
REMC algorithm is properly implemented, and that it 
might give thermalized results for N — 100, n r = 14 up 
the pressure /3P = 38. In this section we study more care- 
fully the outcome of this study, starting with the equation 
of state Z(cp). 

Using the PDFs shown in Fig. [3] it is easy to deduce 
the average volume fraction for each pressure, and thus 
to obtain Z(ip). The results are shown in Fig. [3}a for 
three different systems sizes N = 60, 80 and 100. A 
comparison of the three system sizes shows that finite size 
effects appear to be very small for the equation of state as 
the data obtained with different system sizes practically 
coincide. Nevertheless, larger system sizes produce a very 
small but apparently systematic decrease of the volume 
fraction at a given pressure, for all ip. As a further check, 
we report the results of an independent study using a 
standard Monte-Carlo approach which used N = 1000, 
but covers a smaller range of pressures [44| ■ Up to ip « 
0.595 where both data sets can be compared, the data 
agree very well, confirming the validity of our algorithm, 
at least up to this density. 

A first quantitative result from our study stems from 
data at volume fractions larger than the ones studied 
in Ref. which were all accurately described using 
the Boublik-Mansoori-Carnahan-Starling-Leland (BM- 
CSL) equation of state [56], [5?j], which is the extension 
for mixtures of the Carnahan-Starling equation of state. 
The data in Fig. [5]-a follow the BMCSL equation up to 
ip « 0.59, but clearly deviate from it at larger volume 
fractions, the deviations becoming very large at large ip 
where BMCSL clearly underestimates the pressure. A 
similar deviation was recently reported from molecular 
dynamics simulations of a hard sphere system with con- 
tinuous polydispersity [58[. 

As a better description of the data at large ip we used 
the fitting formula suggested from free volume consider- 
ations [59j |. 



Z 



d'pc 



(10) 



where d! and ip c are free fitting parameters. Although 
the prefactor d! should in principle be constrained within 
free volume theory to be equal to the spatial dimension 
d, we find that its value must be adjusted to describe 
our data. In Fig. [5}a we show the best fit to the data 
ip > 0.61 to Eq. (TTU]) as a full line, using d' = 2.82 and 
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FIG. 5: a- Equation of state Z(ip) for different system sizes, 
N = 60 (open squares), N = 80 (circles), and N = 100 
(diamonds). Light bullets are N = 1000 data taken from 
Ref. [HI]. The black solid line is the BMCSL equation of 
state, the dashed line an empirical polynomial form, and the 
light line a simple pole divergence, Eq. (|10p . b- Isothermal 
compressibility obtained from density fluctuations (symbols) 
or by derivative of the fits shown in panel a. The vertical 
dotted line is at <p v ft, where no thermodynamic signature of 
a glass transition is found. 



ip c = 0.669. This value of tp c should be compared to the 
lower bound for the diverging pressure discussed above 
in Sec. |n] (see Table HJ which was tp iovi = 0.662 go|. The 
large difference between the two values is a direct sign 
that the REMC algorithm has been able to thermalize 
the system much more efficiently. Note also that the fit 
in Eq. (jTDJ) only works at large volume fractions, while at 
low tp values it clearly deviates from the simulation data. 

Finally to account for the crossover region tp « 0.58 — 
0.61 between the BMCSL and free volume fits, we use an 
empirical high order polynomial fit, shown with a dashed 
line. We give no particular emphasis on a physical inter- 
pretation of this fit, which we simply use as a fitting tool 
to obtain the numerical derivative of the pressure, and 
thus the compressibility, in this intermediate Range of 
volume fractions. 

A second important result of our study is obtained 
by considering the vertical line which corresponds to the 
volume fraction (p v ft = 0.615. While an extrapolation of 
the relaxation time divergence using Eq. ([2} indicates the 
possibility of a glass transition occurring at ip v ft, there is 
no corresponding thermodynamic signature in the equa- 
tion of state, in particular no sign that a kink develops 
as the system size increases, at least for the modest N 
values we have been able to study, see Sec. I VIII 

These results are confirmed in Fig.[5]-b which shows the 
evolution of the isothermal compressibility as a function 
of if for the different system sizes. These data are directly 
obtained from the spontaneous density fluctuations, i. e., 
from x = N({p 2 ) — (p) 2 )/(p) 2 , and they directly confirm 
the absence of any jump in the compressibility over this 



range of volume fractions and system sizes, in particular 
near p v f t . 

We also show in this figure the compressibility values 
as obtained from 5p/S(Zp), using the three fits described 
above, namely using the BMCSL equation of state at 
low if (black line), the polynomial fit at intermediate if 
(dashed line) and the free volume fit at high if (light line) . 
There is excellent agreement between both sets of data 
showing that Eq. © is well satified at over the entire 
range of volume fractions. 

The compressibility data simply amplify the results 
obtained for Z(tp). In particular, the good agreement 
between the BMCSL equation of state and the numer- 
ical data is very good up to if ~ 0.56, but deviations 
in fact already appear at moderate volume fractions 
if w 0.57 — 0.59, that are not obvious from the pressure 
itself (see Fig.[S]-a). Similarly, the free volume description 
of the data is only adequate above if — 0.61. These two 
limits make evident the existence of a crossover regime 
if ~ 0.56 — 0.61 where neither approaches work, and the 
only description we have is an empirical polynomial func- 
tion, which, interestingly, shows two changes of the con- 
cavity but no jump. 

To sum up, our simulation data at low and interme- 
diate densities agree with the BMCSL equation of state, 
while at large densities they are much better described 
by a simple divergence at if c = 0.669. The data show no 
jump of x in the studied range of if values, which encom- 
passes both ip mc t and ip v {t fitted dynamic singularities. 



VI. INCREASING THE PRESSURE FURTHER: 
NONEQUILIBRIUM EFFECTS 

In the previous section we found the unexpected result 
that, for modest system sizes, equilibrium data could be 
produced even beyond the fitted location of the VFT 
singularity. In this section we ask whether it is possible 
to go to even larger volume fractions and cross fo — 
0.635, the putative location of the thermodynamic glass 
transition estimated in Ref. [44j . 

To start answering this question we run a simulation 
with N = 60 and a larger number of replicas, n r = 18 
setting the maximum pressure to /3P = 100 and a ge- 
ometric factor of 0.84. As before, we discard the first 
2 x 10 13 trials, and use 2 x 10 13 trial moves to perform 
measurements . 

The results for Z(ip) and x(<P) are shown in Fig. [6j 
while the pressure, contact values of the radial distribu- 
tion functions, and number of neighbors are shown in 
Fig. [71 The data shown in Fig. [6] are consistent with 
those found previously with n r = 14. There is not only a 
good agreement with the data obtained for n r = 14 and 
N = 100, but also with the free volume extrapolation to- 
wards larger tp. Additionally, for ip < 0.63, the measured 
pressure shown in Fig. [TJ-a matches the imposed pressure 
and all structural quantities display a smooth evolution 
with tp, see Fig.[7Jb. Unfortunately, this smooth behavior 
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FIG. 6: Same as Fig. \5\ comparing runs 
14, maximum pressure PP = 38, filled 
n r — 18, maximum pressure pP = 100, 
lines correspond to the BMCSL equation 
free volume fit Eq. (|10l) . the vertical line 
While both data sets coincide below ip 
data at large pressures have not reached 



with (N = 100, n T = 
circles) to (N = 60, 
open circles). Black 
of state, light lines to 
denotes <p v f t = 0.615. 
sa 0.63, the n r = 18 
thermal equilibrium. 




FIG. 7: Same as Fig.[T]for the run with N = 60 and n r = 18. 
The data scatter at large volume fraction, ip > 0.63 indicates 
nonergodic effects. 



is lost for ip > 0.63, see Fig. suggesting that inadequate 
sampling is performed. 

This conclusion is further supported by the data shown 
in Fig. [5] which shows the path in the pressure space 
of three chosen replicas. Despite an acceptance rate for 
replica exchanges being close to 10%, the replicas clearly 
do not sample all thermodynamic states with equal prob- 
ability. In particular, it is clear that the averages at large 
pressure are performed over a very limited number of 
independent configurations, suggesting that an ergodic 
sampling of the phase space is not achieved. Note that 
the third replica, before getting arrested at large pres- 
sure near the end of the run, smoothly travels among the 




5e+12 le+13 1.5e+13 2e+13 2.5e+13 3e+13 
trials 

FIG. 8: Same as Fig.Hfor the run with N = 60 and n r = 18. 
Thermodynamic state '1' corresponds to the highest pressure 
and '18' to the lowest. Contrary to Fig. [2]here the replica do 
not appropriately visit all thermodynamic states in an ergodic 
manner. 



highest 14 pressures, which probably explains why the 
two data sets with n r = 14 and n r = 18 produce con- 
sistent results below (p = 0.63, despite the fact that the 
latter run is clearly not producing fully thermalized data. 

Therefore, we conclude that much longer simulations 
would be needed to reach thermal equilibrium above 
ip ~ 0.63, presumably with a larger number of repli- 
cas to allow a more extensive sampling of the configura- 
tion space. Unfortunately, this implies that despite our 
numerical effort we cannot discuss the possibility raised 
in Ref. [44| that a thermodynamic glass transition takes 
place near ipo = 0.635 in the present binary hard sphere 
mixture. 



VII. DISCUSSION 

In this work, we have demonstrated that the replica ex- 
change Monte-Carlo method recently adapted to improve 
the sampling of hard sphere systems is a useful new tool 
to investigate the thermodynamic behaviour of the disor- 
dered fluid state in a binary mixture of hard spheres up to 
very large volume fractions. We found that reproducible, 
thermalized results could be obtained up to ip ~ 0.63 at 
least for moderate system sizes, N < 100. This vol- 
ume fraction is beyond two important 'critical' packing 
fractions defined dynamically, namely the mode-coupling 
transition (p mc t — 0.592 and the divergence extrapolated 
using a Vogel-Fulcher-Tamman expression, ip v f t — 0.615. 
Following the equation of state for the pressure Z(<p) and 
the isothermal compressibility xif) we found no signa- 
ture of a thermodynamic glass transition up to ip = 0.63 
for the system sizes we use. Additionally, we have pushed 
the lower bound for the location of the divergence of the 
pressure of the fluid branch up to <p — 0.669, much above 
the previous determination ip = 0.662. 

For computational reasons our study was limited both 
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in the range of system sizes and of volume fractions for 
which thermal equilibrium could be reached. Thus, our 
results leave open the existence of (at least) three differ- 
ent scenarios for the behaviour of the fluid of hard spheres 
at large volume fractions in the thermodynamic limit. 

In a first scenario, we assume that finite size effects arc 
small and that our data at large pressure above ip = 0.63 
are nevertheless indicative that no change of behaviour 
is to be expected for Z and x even at larger volume frac- 
tions, such that Eq. (|TU1) will continue to hold up to some 
ip c . In this view, ip c would represent the end point of the 
fluid branch where the equilibrium pressure of the fluid 
would diverge, while the region (p = 0.56 — 0.61 rep- 
resents a crossover from the BMCSL equation of state 
to a free volume-like divergence. To prove or disprove 
this scenario is nearly impossible, as one should establish 
that no thermodynamic singularity occurs up to (p c in 
the thermodynamic limit. It is also natural to expect, 
in this perspective, that the equilibrium relaxation time 
of the fluid would also diverge at u> c . This was termed 
the "jamming" scenario in Ref. [44j because it is the di- 
verging pressure that controls the divergence of the vis- 
cosity. Note that our results imply that this divergence 
will in any case occur above (p = 0.669, which is much 
above the location of the jamming transition ('point J') 
at tpj = 0.648 obtained using purely athermal packing 
preparation protocols [60(. Thus, even in the absence of 
a thermodynamic glass transition, point J does not con- 
trol the glass transition of hard spheres. 

A second scenario could be that finite size effects are 
severe, that our checks with different system sizes are in- 
sufficient, and that N — 100 is still very far away from 
the thermodynamic limit even in the crossover regime 
ip = 0.58—0.62. In that case, the crossover region we have 
described could potentially become sharper, in the ther- 



modynamic limit, yielding a discontinuity of the pressure 
and a jump of the compressibility. This scenario is poten- 
tially simpler to study numerically based on our work, as 
one should attempt to increase further the range of sys- 
tem sizes studied while maintaining thermal equilibrium 
in the crossover regime, an objective that does appear 
numerically realistic. 

In a third scenario, our conclusions would continue to 
hold in the thermodynamic limit up to ip = 0.63, con- 
firming in particular that nothing special happens near 
<Pvft — 0.615. However, a thermodynamic transition 
could still take place at lar ger density, as suggested for 
instance in Refs. [HI, 0, 146U48! where a dynamic singu- 
larity was located near ipo = 0.635. Although we found 
no thermodynamic signature of po in Sec. IVII we also no- 
ticed that our data at these large volume fraction were 
not thermalized leaving open the possibility that a pres- 
sure discontinuity exists at equilibrium. Exploring this 
third scenario would be quite demanding, as one would 
need to cross po at thermal equilibrium for larger sys- 
tems. 

To conclude, it should come as no surprise that pro- 
viding solid conclusions regarding the existence of a ther- 
modynamic liquid-glass transition in the thermodynamic 
limit is a difficult numerical task. However, we have pro- 
vided evidence that replica exchange Monte-Carlo sim- 
ulations can be used to study this issue in hard sphere 
systems, and we have suggested that drawing some firm 
conclusions is perhaps not completely out of reach. In 
particular, it would be interesting to study larger system 
sizes, together with a larger number of replicas to main- 
tain the acceptance rates for replica exchanges at an ac- 
ceptable level. This implies using larger computational 
resources. 
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